/*************************************************************************
Purpose: Loads a FES ascii data file, modifies it in some manner, 
         and writes it to disk using another name.
**************************************************************************
Copyright 2010 Bryan Killett

This file is part of GAIA.

GAIA is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

GAIA is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GAIA. If not, see <http://www.gnu.org/licenses/>.
*************************************************************************/
//Header file "definitions.hpp" declares many functions, constants, datatypes.
//In order to use the functions declared in definitions.hpp, must link
//with object code that includes functions.cpp.
#include "definitions.hpp"

int main(){ 
  // 0 - zero values on land
  // 1 - convert from FES to latlon sin-cos format PROBS AT LON0? 
  // 2 - convert from latlon sin-cos format back to FES. PROBS AT LON0? 
  //10 - Calculate files minus prefix+files in sin-cos format, then report amp/phase diff
  //     2010-07-06 Update: Do this by setting mask values to 0.0 in either file FIRST, to deal with
  //                        situations when FES is masked on land but the other file isn't. Otherwise
  //                        NOTHING is subtracted on land, which is not the intent!
  //11 - Calculate files plus prefix+files in sin-cos format, then report amp/phase sum
  //21 - decimate FES maps in both lat, lon directions by param1.
  //1001 - Replace FES with:
  //       500cm amp, 90 phase in "width" square north of slat, slon.
  //       500cm amp,270 phase in "width" square south of slat, slon.
  //       0cm amp, 0 phase everywhere else.
  //1002 - Replace FES with:
  //       500cm amp, 90 phase in "width" square centered on slat, slon.
  //       0cm amp, 0 phase everywhere else.
  //1003 - Replace FES with:
  //       100cm amp, 0 phase everywhere.
  int choice = 2;

  string prefix = "new/modified-";
  long long param1 = 8;
//      61.0764      293.5714     -2.43325045e+01  -622.96
  double slat=61.0764,slon=293.5714,width=130.0;//Deg,deg,km, used in choice 1001.
  int mode=1;//Used for subtract_fes

  long long i,j,l;
  model_s model;
  latlon_s rawfes, rawfes2;
  model.options.basefolder = "data/FES2004/data/";
  model.options.type=1;
  string basefolder(gaiafolder());
  basefolder.append(model.options.basefolder);

  define_model_titles_and_periods(model);
  
  if(choice==0){
    for(l=0;l<(int)model.options.filenames.size();l++){
      rawfes.filename = basefolder;
      rawfes.filename.append(model.options.filenames.at(l));
      
      rawfes.discard_redundant_values = 1;
      read_rawfes(rawfes);
      
      for(j=0;j<(long long)rawfes.lat.size();j++){
        for(i=0;i<(long long)rawfes.lon.size();i++){
          //Zero both amp and phase if either amp or phase are the mask value, just to be sure
          if(rawfes.amplitudes.at(j).at(i)==rawfes.mask or rawfes.phases.at(j).at(i)==rawfes.mask){
            rawfes.amplitudes.at(j).at(i) = 0.0;
            rawfes.phases.at(j).at(i) = 0.0;
          }
        }
      }

      rawfes.filename = basefolder;
      rawfes.filename.append("new/");
      rawfes.filename.append(model.options.filenames.at(l));

      write_rawfes(rawfes);
    }
  }
  else if(choice==1){
    results_s results;
    grid_s temp_grid;
    plot_options_s plot_options;
    results.options.output_choice = 5;
    plot_options.outputfolder = basefolder;
    plot_options.outputfolder.append("new/");
    string base = "fes_latlon";
    convert_rawfes_to_latlon(results);
    convert_results_from_ampphase_to_sincos(results);
    write_output(base, results, temp_grid, plot_options, 1);
  }
  else if(choice==2){
    results_s results;
    grid_s temp_grid;
    plot_options_s plot_options;
    results.options.output_choice = 5;
    plot_options.inputfolder = model.options.basefolder;
    plot_options.inputfolder.append("new/");
    load_output(results, plot_options, temp_grid);
    convert_latlon_to_rawfes(results);
  }
  else if(choice==10 or choice == 11){
    if(choice==10) mode = 1;
    else if(choice==11) mode = 0;
    for(l=0;l<(int)model.options.filenames.size();l++){
      //Read original FES file.
      rawfes.filename = basefolder;
      rawfes.filename.append(model.options.filenames.at(l));
      rawfes.discard_redundant_values = 1;
      read_rawfes(rawfes);
      
      //Read the file with the prefix in front- hopefully it exists.
      rawfes2.filename = basefolder;
      rawfes2.filename.append(prefix);
      rawfes2.filename.append(model.options.filenames.at(l));
      rawfes2.discard_redundant_values = 1;
      read_rawfes(rawfes2);
    
      //Calculates rawfes-rawfes2 using sin/cos, stores amps/phases in rawfes.
      //When mode==1, this subtracts. When mode==0, it adds.
      subtract_fes(rawfes,rawfes2,mode);

      //Write as a standard ascii FES 2004 file.
      rawfes.filename = basefolder;
      rawfes.filename.append("new/diff-");
      rawfes.filename.append(model.options.filenames.at(l));
      write_rawfes(rawfes);
    }
  }
  else if(choice==21){
    for(l=0;l<(int)model.options.filenames.size();l++){
      rawfes.filename = basefolder;
      rawfes.filename.append(model.options.filenames.at(l));
      
      rawfes.discard_redundant_values = 1;
      read_rawfes(rawfes);

      decimate_fes(rawfes,rawfes2,param1);

      rawfes2.filename = basefolder;
      rawfes2.filename.append("new/");
      rawfes2.filename.append(model.options.filenames.at(l));

      write_rawfes(rawfes2);
    }
  }
  else if(choice==1001){
    for(l=0;l<(int)model.options.filenames.size();l++){
      rawfes.filename = basefolder;
      rawfes.filename.append(model.options.filenames.at(l));
      
      rawfes.discard_redundant_values = 1;
      read_rawfes(rawfes);

      double mid_lat = slat;
      double mid_lon = slon;
      double width_lat = width*km2deg;
      double width_lon = width_lat/cos(deg2rad*mid_lat);
      double min_lat = mid_lat-width_lat;
      double max_lat = mid_lat+width_lat;
      double min_lon = mid_lon-width_lon/2;
      double max_lon = mid_lon+width_lon/2;
      
      for(j=0;j<(long long)rawfes.lat.size();j++){
        for(i=0;i<(long long)rawfes.lon.size();i++){
          if((rawfes.lat.at(j) > mid_lat and rawfes.lat.at(j) <= max_lat) and (rawfes.lon.at(i) > min_lon and rawfes.lon.at(i) <= max_lon)){
            rawfes.amplitudes.at(j).at(i) = 500.0;
            rawfes.phases.at(j).at(i) = 90.0;
          }
          else if((rawfes.lat.at(j) > min_lat and rawfes.lat.at(j) <= mid_lat) and (rawfes.lon.at(i) > min_lon and rawfes.lon.at(i) <= max_lon)){
            rawfes.amplitudes.at(j).at(i) = 500.0;
            rawfes.phases.at(j).at(i) = 270.0;
          }
          else{
            rawfes.amplitudes.at(j).at(i) = 0.0;
            rawfes.phases.at(j).at(i) = 0.0;
          }
        }
      }

      rawfes.filename = basefolder;
      rawfes.filename.append("new/");
      rawfes.filename.append(model.options.filenames.at(l));

      write_rawfes(rawfes);
    }
  }
  else if(choice==1002){
    for(l=0;l<(int)model.options.filenames.size();l++){
      rawfes.filename = basefolder;
      rawfes.filename.append(model.options.filenames.at(l));
      
      rawfes.discard_redundant_values = 1;
      read_rawfes(rawfes);

      double width_lat = width*km2deg;
      double width_lon = width_lat/cos(deg2rad*slat);
      double min_lat = slat-width_lat/2;
      double max_lat = slat+width_lat/2;
      double min_lon = slon-width_lon/2;
      double max_lon = slon+width_lon/2;
      
      for(j=0;j<(long long)rawfes.lat.size();j++){
        for(i=0;i<(long long)rawfes.lon.size();i++){
          if((rawfes.lat.at(j) > min_lat and rawfes.lat.at(j) <= max_lat) and (rawfes.lon.at(i) > min_lon and rawfes.lon.at(i) <= max_lon)){
            rawfes.amplitudes.at(j).at(i) = 500.0;
            rawfes.phases.at(j).at(i) = 90.0;
          }
          else{
            rawfes.amplitudes.at(j).at(i) = 0.0;
            rawfes.phases.at(j).at(i) = 0.0;
          }
        }
      }

      rawfes.filename = basefolder;
      rawfes.filename.append("new/");
      rawfes.filename.append(model.options.filenames.at(l));

      write_rawfes(rawfes);
    }
  }
  else if(choice==1003){
    for(l=0;l<(int)model.options.filenames.size();l++){
      rawfes.filename = basefolder;
      rawfes.filename.append(model.options.filenames.at(l));
      
      rawfes.discard_redundant_values = 1;
      read_rawfes(rawfes);

      for(j=0;j<(long long)rawfes.lat.size();j++){
        for(i=0;i<(long long)rawfes.lon.size();i++){
          rawfes.amplitudes.at(j).at(i) = 100.0;
          rawfes.phases.at(j).at(i) = 0.0;
        }
      }

      rawfes.filename = basefolder;
      rawfes.filename.append("new/");
      rawfes.filename.append(model.options.filenames.at(l));

      write_rawfes(rawfes);
    }
  }
  else cout<<"Choice "<<choice<<" wasn't recognized!"<<endl;
  return 0;
}
